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ABSTRACT 

In the recent years, the "Nice" model of solar system formation has attained an unprecedented level 
of success in reproducing much of the observed orbital architecture of the solar system by evolving 
the planets to their current locations from a more compact configuration. Within the context of this 
model, the formation of the classical Kuiper belt requires a phase during which the ice giants have a 
high eccentricity. An outstanding question of this model is the initial configuration from which the 
Solar System started out. Recent work has shown that multi-resonant initial conditions can serve as 
good candidates, as they naturally prevent vigorous type-II migration. In this paper, we use analytical 
arguments, as well as self-consistent numerical N-body simulations to identify fully-resonant initial 
conditions, whose dynamical evolution is characterized by an eccentric phase of the ice-giants, as well 
as planetary scattering. We find a total of eight such initial conditions. Four of these primordial states 
are compatible with the canonical "Nice" model, while the others imply slightly different evolutions. 
The results presented here should prove useful in further development of a comprehensive model for 
solar system formation. 

Subject headings: celestial mechanics — planets and satellites: formation — methods: analytical — 
methods: numerical 



1. INTRODUCTION 

The question of how the solar system formed dates 
back centuries. The last decade, however, has seen a 
considerable amount of progress made on this issue. No- 
tably, the development of the Nice model (Tsiaganis et al. 
2005) has proven to be a milestone. The scenario foretold 
by the Nice model is as follows: the giant planets form 
in a compact configuration, and driven by planetesimal 
scattering (Fernandez & Ip 1984, Malhotra 1995) begin 
migrating divergently. Eventually, Jupiter and Saturn 
cross their mutual 2:1 mean-motion resonance (MMR), 
which results in an acquisition of eccentricities for both 
planets. Subsequently, the whole outer solar system un- 
dergoes a brief period of dynamical instability, during 
which Uranus and Neptune are scattered to their cur- 
rent orbits. 

There are a few aspects to the success of the Nice 
model. First and foremost, it has been able to replicate 
the architecture of the secular dynamics of the outer so- 
lar system (Morbidelli et al. 2009a). Second, it provides 
a semi-quantitative description of the formation of the 
Kuiper belt (Levison et al. 2008). Third, the inward flux 
of planetecimals during the phase of dynamical instabil- 
ity allows for chaotic capture of Jupiter's and Neptune's 
Trojan populations (Morbidelli et al. 2005, Nesvorniy et 
al. 2007). Finally, if the resonance crossing between 
Jupiter and Saturn is timed appropriately, the global 
mayhem provides a natural trigger for Late Heavy Bom- 
bardement (LHB) (Gomes et al. 2005). There are other 
observational constraints that should be reproduced in a 
model for the solar system's formation, such as the dy- 
namical structure of the inner solar system, and consid- 
erable progress has already been made in this direction 
(Brasser et al. 2009, Morbidelli et al. 2009b). At the 
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same time, it is also crucial to explore the unobservable 
aspect of the model, namely the initial conditions. This 
is the purpose of our study. 

In this paper, we consider various multi-resonant con- 
figurations as possible initial conditions for the Nice 
model. We investigate the early stages of dynamical evo- 
lution for a large number of candidate systems and show 
that only eight configurations appear to be consistent 
with the formation of the Kuiper belt in the framework 
of the Nice model. It is noteworthy that running simula- 
tions of the Nice model to completion is computationally 
expensive. Consequently, as we seek to examine a large 
array of initial conditions, we are forced to utilize early 
dynamical events, namely planetary scattering, as prox- 
ies for successful formation of the solar system. In this 
manner, we limit the duration of each simulation to only 
a few tens of millions of years. Within the context of our 
integrations, this is long enough for the system to pass 
through the epoch of dynamical instability, but not long 
enough to scatter away all of the planetesimals which end 
up on long-term unstable orbits. 

The plan of the paper is as follows: in section 2, we 
explain how multi-resonant configurations prevent type- 
II migration, and our approach to their assembly. In 
section 3, we discuss the evolution scenarios of the con- 
sidered systems, estimate the amplitudes of eccentricity 
jumps in relevant cases, and present the results of N- 
body simulations. We discuss our results and conclude 
in section 4. 

2. MULTI-RESONANT CONFIGURATIONS 

One of the important differences between the solar sys- 
tem and the majority of the detected aggregate of extra- 
solar planetary systems is the lack of a close-in giant 
planet. This difference suggests that while it is com- 
mon for planets to migrate to small orbital radii, some 
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mechanism was at play in the early solar system which 
prevented vigorous orbital decay. One such mechanism, 
which is both efficient and reasonable, is resonant cap- 
ture (Masset k Snellgrove 2001). 

When a newly formed gaseous planet reaches a critical 
mass of ~ lMj, it opens a gap in the proto-planeatary 
disk. Incidentally, the planet continues to interact with 
the disk via various resonances. Summed together, the 
resonant torques from a given side of the disk, some- 
what counter-intuitively, push the planet away from that 
side. As a result, the planet positions itself at a point 
in the gap where all torques cancel, and moves inward 
together with the disk on the viscous time-scale, in a 
process termed Type-II migration (Morbidelli & Crida 
2007). Simultaneously if another planet, whose semi- 
major axis is larger, is migrating inwards faster, it will 
eventually encounter a mean-motion resonance. Under 
a large spectrum of circumstances, converging orbits can 
lead to capture into a mean-motion resonance, ensuring 
that the two planets' period ratio remains constant for 
extended periods of time. In fact, for slow enough mi- 
gration rates and low enough eccentricities, resonant cap- 
ture is certain (Peale 1986). When this happens, the gas 
between the two planets drains as the gaps overlap. Con- 
sequently, the torque balance on the resonant pair results 
from gas interior to the inner planet and that exterior to 
the outer planet. This leads to a drastic reduction of the 
migration rate (Lee & Peale 2002). Furthermore, if the 
inner planet is more massive than the outer planet, as is 
the case with Jupiter and Saturn, the migration of the 
resonant pair can be halted altogether or even reversed 
(Morbidelli & Crida 2007). 

Numerical studies of Jupiter and Saturn submerged 
in a gaseous proto-planetary disk suggest precisely the 
above scenario. Indeed, Saturn's migration is consider- 
ably faster than Jupiter's because of its lower mass and 
inability to fully open a clean gap. The pioneering re- 
sults of Masset & Snellgrove (2001) showed that locking 
Jupiter and Saturn in the 3:2 MMR can effectively halt 
the pair's migration. The somewhat more precise nu- 
merical experiments of Morbidelli & Crida (2007) con- 
firmed this and also showed that capture into 2:1 and 
5:3^ MMR's are viable outcomes, depending on where 
Saturn forms relative to Jupiter. The work of Pierens & 
Nelson (2008) however suggests that while capture into 
2:1 and 5:3 MMR's is certainly possible, in a number of 
cases Saturn eventually breaks away and continues its 
inward migration until it is captured in the 3:2 MMR. 
Collectively, the above mentioned results suggest that 3:2 
MMR is indeed a likely initial configuration of Jupiter & 
Saturn, although there is not enough evidence to deci- 
sively rule out the 2:1 MMR or the 5:3 MMR as initial 
conditions. Consequently, for the sake of completeness, 
we consider all three of these resonances as possible start- 
ing configuratoins for Jupiter and Saturn. 

By extension of the above scenario, the ice giants, 
which are believed to form after Jupiter and Saturn, 
behave in a qualitatively similar way. Namely, as they 
migrate from the outer disk inwards, they too become 

1 It must be noted that capture into 5:3 MMR is less probable, 
since it is a second order resonance. Furthermore, even if Jupiter 
and Saturn are captured, subsequent motion can be unstable (Mor- 
bidelli & Crida 2007). 



trapped in MMR's. Consequently, at the epoch of the 
disappearance of the gas, we are left with a multi- 
resonant system, in which each planet is in resonance 
with its neighbors. Morbidelli et al. (2007) performed 
hydrodynamical simulations of this process with the pur- 
pose of identifying such configurations that are long-term 
stable. Considering only systems where Jupiter and Sat- 
urn are locked in a 3:2 MMR, they were able to find two 
stable fully-resonant states. 

Our approach to assembling multi-resonant systems 
follows that of Lee & Peale (2002). In addition to Newto- 
nian N-body interactions, each planet is subject to semi- 
major axis decay 



and eccentricity damping 



where a is semi-major axis, e is the eccentricity, and 
K is an adjustable migration frequency. In our simu- 
lations, we keep K the same for all planets, ensuring al- 
ways convergent migration. In accord with Lee & Peale 
(2002), a Bulirsch-Stoer integration method (Press et al. 
1992) was used. In contrast with the full hydrodynami- 
cal simulations, this method is simpler and computation- 
ally cheaper, allowing us to sample a large array of sys- 
tems. Additionally, given the problem's straight-forward 
nature, it is unlikely that a configuration found using this 
approach cannot be obtained using other methods. 

In the context of these simulations, all four giant plan- 
ets were introduced simultaneously on planar circular or- 
bits, slightly outside of their desired resonant locations, 
with the more massive ice giant on the outermost orbit. 
If capture into the desired resonances did not occur, we 
varied K. As pointed out in Morbidelli et al. (2007), the 
sequence in which planets get captured may be impor- 
tant, since changing the order can change the librating 
resonant angles. To avoid this degeneracy, we always set 
the initial orbits such that Saturn would get captured 
first and Neptune last. To avoid confusion, we shall al- 
ways refer to the outermost planet as Neptune, although 
in the Nice-model, the orbits of the ice giants may switch 
places. After each desired configuration was achieved, 
the 5-body system was subjected to a lOOMyr dynami- 
cal stability test, using the mercury6 software (Chambers 
1999). Note that lOOMyr is highly conservative, given 
that our full dynamical evolution simulations only last 
~30Myr. However, these stability integrations show that 
the same initial conditions are also applicable for scenar- 
ios where the global instability occurs somewhat later 
than what is considered in this work. Table (1) lists all 
stable multi-resonant configurations that we generated. 

3. DYNAMICAL EVOLUTION 

Having found a large array of stable multi-resonant 
systems, we now need to determine which of these con- 
figurations can resemble the current state of the solar 
system, having dynamically evolved. There are two con- 
straints of interest here. First is the structure of the 
secular dynamics of the giant planets. Recently, Mor- 
bidelli et al. (2009a) showed that a smooth migration 
scenario, such as the one envisioned by Malhotra (1995), 
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TABLE 1 
Multi-resonant Initial Conditions 
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5:3 


2:1 


4:3 


5:3 


2:1 


5:4 


5:3 


2:1 


6:5 
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4:3 
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3:2 


3:2 


2:1 


3:2 
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2:1 


4:3 


3:2 


2:1 


4:3 


4:3 



All of the stable multi-resonant initial conditions considered in 
this study. The bold lines represent the configurations that proved 
to be compatible with a Nice model-like evolution. 

is incompatible with the observed eccentricities and in- 
clinations of the giant planets. A mean motion resonance 
crossing event by itself is also insufficient because it does 
not excite the inclinations to a necessary degree or repro- 
duce the amplitudes of the g 5 and g 6 secular eigenmodes 
correctly. To create the current eccentricities, inclina- 
tions, and eigenmode amplitudes, encounters must have 
happened between an ice giant and a gas giant. Further- 
more, if the instability took place after the inner solar 
system was already intact, encounters must have taken 
place between an ice giant and both gas giants to cause 
Jupiter to "jump", in order to prevent slow secular reso- 
nance sweeping of the inner solar system (Brasser et al. 
2009). 

The second criterion of interest is the formation of the 
Kuiper belt, particularly the classical region. The trans- 
port mechanism of planetesimals to this region, proposed 
by Levison et al. (2008) relies on overlapping mean mo- 
tion resonances. Namely, when Neptune's eccentricity 
exceeds ~ 0.2, its exterior MMR's widen enough to over- 
lap, and motion of all particles in the region becomes 
highly chaotic. This allows for planetesimals to execute 
a random-walk and invade the classical Kuiper belt re- 
gion. With time, as Neptune's eccentricity decays due to 
dynamical friction (Stewart & Withrill 1988), the reso- 
nances become narrower, and the particles occupying the 
classical region, no longer chaotic, remain on their orbits 
forever. 

Using the two constraints described above as proxies 
for a successful formation scenario, we look for a subset of 
our generated initial conditions that result in evolutions 
which encompass both a scattering event and a transient 
high-eccentricity phase of the outer ice giant. In the Nice 
model, these two constraints are practically always sat- 



isfied simultaneously, since the ice giants tend to switch 
places and the scattering event is intimately tied to the 
high eccentricities. 

Our dynamical evolution simulations include the four 
outer planets and a disk of ~ 3000 equal-mass planetes- 
imals, while the mass of the inner solar system is added 
to the sun. The radial surface density of the planetesi- 
mal swarm was assumed to have a power law structure: 
E oc r~ k , where k £ (1,2). Consequently, the set of 
initial-value problems at hand is controlled by three pa- 
rameters: Jupiter's semi-major axis, which due to reso- 
nant relations controls the semi-major axes of the other 
planets, the planetesimal disk's mass, rridisk, arL d k. As it 
turns out, the actual value of k has little effect on whether 
a given initial condition gives rise to a scattering event 
and an eccentric Neptune. Rather, it controls how fast 
the planets grind through the disk. As a result, we allow 
it to float randomly from simulation to simulation. Also, 
an advantage of multi-resonant initial conditions lies in 
that with enough simulations, it is possible to determine 
a unique combination of (aj,mdisk), since Jupiter and 
Saturn must scatter enough particles to arrive to their 
current 5:2 commensurability. 

In our simulations, the inner edge of the planetesimal 
disk was placed <~ 1AU outside of Neptune's orbit. This 
forces migration, driven by planetesimal scattering, to 
begin shortly after the start of the simulation. In other 
words, we do not attempt to time the onset of instability 
with LHB, as was done in Gomes et al. (2005). There 
is another implication of a relatively close inner edge. If 
the planetesimal swarm is nearby when the instability be- 
gins, the planetary orbits penetrate deeply into the disk, 
and the resulting dynamical friction plays a stabilizing 
effect (Levison, personal communication). Intuitively, it 
makes sense to place the inner edge of the disk where the 
dynamical lifetime of planetesimals equals the lifetime of 
the gaseous nebula. While a 1 AU separation is approxi- 
mately correct (Gomes et al. 2005), in the future it may 
be a worthwhile exercise to determine this boundary pre- 
cisely for each multi-resonant initial condition. The outer 
edge of the disk was placed at 30 AU with the purpose 
of eventually halting Neptune's migration (Tsiaganis et 
al. 2005, Levison et al. 2008). 

In order to not force the scattering-driven migration 
fictitiously, but still keep the computational cost down, 
we ignore self-gravity of the planetesimal disk. The re- 
sulting Hamiltonian takes the form: 

i=l 1 i=l j=i+l 1 n 

where p is momentum, x is position, m is mass and G is 
the universal gravitational constant. To further dimin- 
ish the computational cost, after a particle was scattered 
beyond 500AU, it was removed from the simulation. It 
must be noted that planetcsimal-planetesimal interac- 
tions may in reality be important, since they give rise 
to an effective viscosity in the swarm (Levison, personal 
communication). Thus, more careful validation of our 
results should preferentially include these thorny effects. 
We used a hybrid Bulisch-Stoer/Wisdom-Holman algo- 
rithm of the mercury6 software for all integrations. We 
consistently used a time-step of r = 300d, and checked all 
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successful simulations with a smaller r = 60d time-step 
to ensure that the observed instability is not a numerical 
artifact. In all such checks, the evolutions were prac- 
tically indistinguishable, which assures that integrals of 
motion are sufficiently conserved. Finally, our simula- 
tions only cover a few tens of millions of years, since the 
focus here is on distinguishing between initial conditions 
that give rise to gas-giant /ice-giant scattering and ones 
that don't. As a result, the long-term evolution of the 
system after the instability is unexplored. 

3.1. Initial Conditions with Jupiter and Saturn in a 3:2 

MMR 

Let us first consider a family of initial conditions, listed 
in table (1), where Jupiter and Saturn are in a 3:2 MMR. 
This family of initial conditions was previously studied 
in some detail by Morbidelli et al. (2007). Consequently, 
this section's results are partially reproductions. Using 
a hydrodynamical model, Morbidelli et al. (2007) found 
six multi-resonant configurations, two of which they de- 
termined to be long-term stable. These are the configura- 
tions listed in table (1) where both Jupiter & Saturn and 
Saturn & Uranus pairs are in 3:2 MMR's while Uranus & 
Neptune are in either 4:3 or 5:4 MMR's. There are two 
more compact configurations listed in table (1) which 
we determined to be stable, although the counterparts 
of these configurations put together by Morbidelli et al. 
(2007) were unstabl^J These configurations are the two 
where Jupiter & Saturn are in a 3:2 MMR, Saturn & 
Uranus are in a 4:3 MMR, and Uranus & Neptune are 
either in a 3:2 MMR or 4:3 MMR. 

As discussed in Morbidelli et al. (2007), if Jupiter and 
Saturn start out in a 3:2 MMR, the instability is trig- 
gered by their encounter with the 5:3 MMR. While this 
is a second-order resonance, small jumps in Jupiter's and 
Saturn's eccentricities go a long way, especially in highly 
compact configurations. Unfortunately, in this case it is 
difficult to conclusively determine which configurations 
will result in evolutions with scattering events a-priori. 
Thus, we must rely solely on numerical integrations to 
explore the various evolutionary outcomes of these ini- 
tial conditions. 

After an initial run of 20 integrations for each initial 
condition of the family listed in table (1), we ruled out 
the configurations where Saturn & Uranus are in a 2:1 
MMR as well as the configuration where all planet pairs 
are in 3:2 MMR's because all evolutions were character- 
ized by smooth migration. We subjected the remain- 
ing four configurations to 30 additional integrations and 
found that the only configuration which does not result 
in ice-giant /gas-giant scattering is the one where Saturn 
& Uranus are in a 3:2 MMR and Uranus & Neptune are 
in a 4:3 MMR. The evolutions of the remaining initial 
conditions are presented in figures (1) - (3), and their 
final orbital parameters are entered into table (2). 

For the initial condition in which Saturn and Uranus 
are in a 3:2 MMR, 10% of the integrations were successful 
with 57% of them exhibiting close encounters between an 
ice giant and both gas giants. The same fractions for the 
two configurations where Saturn and Uranus are initially 

2 Note that our simplified migration model is computationally 
cheaper than that of Morbidelli et al. (2007), allowing us to run 
more trials to arrive at stable configurations. 
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Fig. 1. — Dynamical evolution of the initial configuration where 
initially Jupiter & Saturn are in a 3:2 MMR, Saturn & Uranus 
are in a 3:2 MMR and Uranus & Neptune are in a 5:4 MMR (as 
labeled). Curves depicting the semi-major axes, perihelion, and 
apohelion of each planet are labeled. The current semi-major axes, 
perihelia and apohelia of the current solar system are plotted as 
grey points for comparison. In this model, the global instability 
is brought forth by Jupiter & Uranus encountering a mutual 5:3 
MMR. 




a 4:3 MMR and Uranus & Neptune are in a 3:2 MMR. 

in a 4:3 MMR are 20% & 30% and 27% & 50% for the 
case where Uranus $ Neptune are in a 3:2 MMR and 4:3 
MMR respectively. 

3.2. Initial Conditions with Jupiter and Saturn in a 5:3 

MMR 

We now move on to the next family of initial condi- 
tions. To begin with, we take the same approach as 
above. Stable multi-resonant configurations of this fam- 
ily are listed as the second set of entries in table (1). We 
simulated the evolutions of these systems with 20 integra- 
tions each. After completion, a clear boundary between 
initial conditions that result in smooth migration and 
those that result in scattering developed. Namely, all se- 
tups where Saturn & Uranus are initially in a 2:1 MMR 
were characterized by smooth evolutions. A similar sce- 
nario describes the fate of initial conditions where Saturn 
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Fig. 3. — Same as fig.l except initially, both Saturn & Uranus 
and Uranus & Neptune are in 4:3 MMR's. 
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Fig. 4. — Dynamical evolution of the initial configuration where 
initially Jupiter &: Saturn are in a 5:3 MMR, Saturn & Uranus 
are in a 3:2 MMR and Uranus & Neptune are in a 5:4 MMR. In 
this model, the global instability is brought forth by Jupiter & 
Saturn encountering a mutual 2:1 MMR, just as in the classical 
Nice model. All else is as in fig.l. 

& Uranus are in a 3:2 MMR while Uranus & Neptune are 
in a 3:2 or a 4:3 MMR. However in the same context, if 
Uranus and Neptune start out in a 5:4 or a 6:5 MMR, 
ice giant/gas giant scattering as well as transient phases 
of high eccentricities are present. Particularly, for the 
configuration where Uranus & Neptune start out in a 
5:4 MMR, 20% of the integrations were successful with 
50% of them exhibiting close encounters between an ice 
giant and both gas giants. For the configuration where 
Uranus & Neptune start out in a 6:5 MMR, also 20% of 
the integrations were successful, but none of the solutions 
exhibited ice giant encounters with both gas giants. 

In the subset of initial conditions where Saturn & 
Uranus are in a 4:3 MMR, the configurations with Uranus 
& Neptune in a 3:2 MMR and a 4:3 MMR can serve as 
good candidates for solar system formation, but the con- 
figuration with Uranus & Neptune in a 5:4 MMR con- 
sistently leads to ejections. 10% of the integrations with 
Uranus & Neptune initially in a 3:2 MMR were success- 
ful, all of them exhibiting scattering with both ice giants. 



In the context of the configuration with Uranus & Nep- 
tune initially in a 4:3 MMR, 15% of the integrations were 
successful, while 33% of them lead to encounters of an 
ice giant with both gas giants. 

Examples of successful evolutions that start from the 
initial conditions described above are presented in figures 
(4) - (7), with final orbital parameters entered into table 
(2). Note that the scattering event in the evolution of 
the configuration where Saturn & Uranus are initially in 
a 3:2 MMR and Uranus & Neptune are in a 5:4 MMR 
(fig. 4) is considerably more violent than that in most 
other examples. This is because majority of close en- 
counters here is between Jupiter and Uranus, while in 
most other simulations, Saturn is responsible for scat- 
tering. It is important to understand that this is not 
a unique feature of the particular initial condition. We 
have observed similar phenomena in simulations of other 
setups as well. 

Note that in a scenario where Jupiter and Saturn start 
out in a 5:3 MMR, the instability is brought on by their 
crossing of the 2:1 MMR, just as in the classical Nice 
model. Much effort has been put into fine-tuning the 
classical Nice model's initial conditions (Tsiaganis et al. 
2005, Morbidelli et al. 2005, Gomes et al. 2005, Levison 
et al. 2008). What matters most, however, are the loca- 
tions of the planets when Jupiter and Saturn are crossing 
the 2:1 MMR. Let us now examine if the classical Nice 
model is compatible with any multi-resonant initial con- 
ditions from this family. 

We begin our calculation by measuring the semi-major 
axes of the four planets at the Jupiter/Saturn 2:1 MMR 
crossing in figure (2a) of Levison et al. (2008). The val- 
ues are listed in the second column of Table (3). Between 
encounters with MMR's, migrations of Jupiter and Sat- 
urn are mostly due to scattering of planetesimals. Mal- 
hotra (1995) showed from conservation of angular mo- 
mentum that this process obeys a relation, which upon 
integration can be written as 

Aloga^/^^- (4) 

TYlplanet 

where / is an empirically determined "efficiency," listed 
in Table (3), and m s is the total scattered mass. When 
applied to Jupiter and Saturn simultaneously, this re- 
lation can be used to "backtrace" the system's migra- 
tion, by roughly estimating the starting semi-major axes 
of Jupiter and the scattered mass m s . Setting Aas — 
(5/3) 2 /V 7 - {2) 2 / 3 a f j and Aa.j = a f j- a l j with a f j = 5.45 
AU yields a total scattered mass of 37m® and a % j = 5.69 
AU. 

Neglecting high-order resonant encounters, we then use 
the calculated scattered mass and apply equation (4) to 
Uranus and Neptune to determine their original posi- 
tions. The back-traced initial conditions are listed in the 
third column of table (3). Incidentally, these initial con- 
ditions are close to a multi-resonant configuration where 
Saturn & Uranus and Uranus & Neptune are both in 4:3 
MMR's. Recall that this initial condition is indeed one of 
the setups that consistently exhibit scattering. However, 
given the similarities in dynamical evolutions among the 
successful initial conditions of this family, at this level of 
accuracy, it is probably safe to say that all four of them 
are compatible with the classical Nice model results. 
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Fig. 5. — Same as fig. 4 except initially, Saturn & Uranus are in 
a 3:2 MMR and Uranus & Neptune are in a 6:5 MMR. 
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Fig. 6. — Same as fig. 4 except initially, Saturn & Uranus are in 
a 4:3 MMR and Uranus & Neptune are in a 3:2 MMR. 
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Fig. 7. — Same as fig. 4 except initially, both Saturn & Uranus 
and Uranus & Neptune are in 4:3 MMR's. 



3.3. Initial Conditions with Jupiter and Saturn in a 2:1 

MMR 

Let us now consider the final family of initial condi- 
tions, listed in table (1), where Jupiter and Saturn are 
initially in a 2:1 MMR. Unlike the scenario of the classi- 
cal Nice model (Tsiaganis et al. 2005), there are no ma- 
jor resonances to cross for Jupiter and Saturn between 
the 2:1 and the 5:2 MMR's. Consequently, a different 
mechanism, involving different resonances, is needed to 
create the instability. Thommes et al. (2008) consid- 
ered the dynamical evolution of a system where Jupiter 
& Saturn are in a 2:1 MMR, Saturn & Uranus are in a 
3:2 MMR, and Uranus & Neptune are in a 4:3 MMR. In 
such a system, the instability is triggered by Uranus and 
Neptune crossing a 7:5 MMR. Due to a weaker, second- 
order nature of this resonance, the eccentricity increase is 
rather small. Incidentally in this particular system, this 
is enough for the ice giants to cross orbits and scatter off 
of each other, but not off of one of the gas giant^] It 
appears that somewhat larger eccentricities are needed. 
Testing each initial condition with a large number of nu- 
merical simulations, as discussed above, is rather time- 
consuming. Consequently, it is worthwhile to quantify 
the amplitudes of eccentricity jumps due to various res- 
onance crossings before-hand if possible. For this set of 
initial conditions, under the assumption of adiabatic mi- 
gration, the eccentricity jumps are deterministic and can 
be estimated analytically (Henrard 1982). 

Following the treatment of (Peale 1986, see also Mur- 
ray & Dermott 1999), we consider the planar internal 
first-order j : (j — 1) resonant Hamiltonian 

G 2 m,r,m 3 G 2 m,r,m' 3 

TT _ O 



2A 2 
"rngTOTO" 

■ r7 sec + AA, 



2A' 2 



7(a/«')y T cos(jA' + (l-j)A + 7 ) 

3C - ^'iLc + A'A' sec 



where A is the mean longitude, 7 



(5) 

vj is the longitude 



of perihelion, A = (m m Q )/(r7i + mQ)y/G(mQ + m)a 

& r = A(l — Vl - e 2 ) are their respective Poincare 
conjugate momenta, and the prime designates the outer 
planet. The secular changes in mean longitude and lon- 
gitude of perihelion are accounted for by the last four 
terms, while f(a/a') arises from the classical expansion 
of the planetary disturbing potential and is a function 
of Laplace coefficients and their derivatives. The expres- 
sions for f(a/a') are presented in Appendix B of (Murray 
& Dermott 1999). Under a series of variable transforma- 
tions (see Peale 1986 for derivation), this Hamiltonian 
can be rewritten to take a simpler form. Let us intro- 
duce the constants a. 8 and e: 



a 



(j-l)n*- jn'* + % 



(J I) 2 



(6) 



(7) 



3 Another scenario present in the integrations of Thommes et 
al. (2008) is the escape of one of the ice giants. In this case, the 
remaining ice giant is left with a high eccentricity, but there are 
only 3 planets left in the system. 
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= « 3 / 2 , 



2 / 

a m 



' f{a/ a ')—r Vm, 

a' m & 



(8) 



where n* is sum of the Keplerian mean motion and the 
secular change in mean longitudes. It is important to 
note that these expressions are not strictly constant, 
since semi-major axis changes. However, (3 is only weakly 
dependent on semi- major axis (Pealc 1976), and in the 
case of e, variations due to the cosine term dominate, so 
the assumption of constant coefficients is sound (Murray 
and Dermott 1999). Relative to the original Hamilto- 
nian, we scale the momentum as: 



2/3 



2/3 



(9) 



The corresponding conjugate angle <j) is simply the cosine 
argument in equation (5), although if e > 0, we also need 
to add 7r to the expression (Murray & Dermott 1999). 
That said, the transformed Hamiltonian takes the form 

(10) 



H res = 5<5> + $ 2 - 2V2$ cos(0). 
This Hamiltonian is parameterized by 



a 



e 2 /3 



1/3 



(11) 



which is a measure of the perturbed object's proximity 
to exact resonance. Finally we note that this Hamil- 
tonian is most easily visualized in terms of polar coordi- 
nates, so we introduce the new mixed canonical variables 
x = \/2i> cos 4> and y = V2<$> sin <f> (Henrard & Lamaitre 
1983). The Hamiltonian now takes the form 



S(x 2 + y 2 ) (x 2 +y 2 ) 



2\2 



- 2x. 



(12) 



Upon application of Hamilton's equations, we see that 
the stationary points of the above Hamiltonian are de- 
scribed by the equation 



+ Sx - 2 = 0. 



(13) 



For resonant encounters aided by divergent migration, 
S < 3 initially. In this case, the existence of a separatrix 
is ensured, and there are three real fixed points, all of 
which lie on the x-axis. Two of these points are always 
negative, and the more negative one is unstable, as it lies 
on the intersection of the inner and the outer branches 
of the separatrix. This is crucial to the estimation of 
eccentricity jumps during resonant encounters. 

If migration is slow enough for 5 to be approximately 
constant over one period of motion, the action, defined 
as 



J 



= <j> <&d(j) = <j> 



xdy 



(14) 



is an adiabatic invariant (Pcale 1986). In other words, it 
is constant except during separatrix crossing. Further- 
more, when the separatrix is far away, the trajectories of 
the circulating orbits in (x, y) space are circles to a good 
approximation. Consequently, we can write J = 2ir§ 
(Murray & Dermott 1999). 

When two planets approach commensurability, a wide 
separatrix is seen as shrinking down on the orbit of the 
perturbed planet in (x, y) space. When the inner branch 



of the separatrix engulfs the planetary orbit, the pro- 
cess of resonance crossing is characterized by the planet 
switching to the separatrix's outer circulating branch. 
The outer branch has a wider radius, thus the increase 
in action. However, during this switch, the perturbed 
planet must necessarily pass through the unstable sta- 
tionary point described above. Consequently, the cal- 
culation is as follows: knowing the action prior to the 
resonant encounter, we can determine the value of 5 at 
the transition using equation (13). Recall however, that 
S also parameterizes the Hamiltonian, and therefore de- 
termines the shape of the separatrix, while the area en- 
gulfed by the outer branch corresponds to the new action 
(see supplemental material of Tsiaganis et al. 2005 for 
an intuitive discussion). It can be shown that the actions 
before and after resonance crossing are related by 



J t + J f = -2ir5. 



(15) 



Thus, the new eccentricity can be easily backed out. 

The above analysis can also be applied to external res- 
onances. In this case, 7 in the cosine argument of the 
Hamiltonian (5) is replaced by 7', and its factor y / 2T T /A 
is replaced by y/2T' / A' , since we are now concerned with 
an e! resonance. Accordingly, we change the scaling fac- 
tors to 

a=(j-l)n*-jn'*+^ ec , (16) 



i' 3 / 2 f(a/a')a'- 



<m 



(17) 



while P remains the same. Note also that any indirect 
terms in the expansion of the disturbing function must be 
accounted for in f(a/a'). Under these transformations, 
Hamiltonian (10) still applies, and so does the subsequent 
analysis (Murray and Dermott 1999). 

The resulting estimates of eccentricity jumps for var- 
ious first-order resonant encounters between Uranus & 
Neptune and Saturn & Uranus are listed in Table (4). 
As can be seen from these calculations, all first-order res- 
onant encounters between Uranus and Neptune produce 
rather small eccentricity jumps. Therefore, we disfavor 
them as good options for triggering instability scenarios 
in which encounters with Saturn take place. Numerical 
integrations performed in the two previous sections arc 
suggestive of this as well. We therefore rule out the con- 
figurations where Saturn and Uranus are in a 2:1 MMR. 
Resonant encounters between Saturn and Uranus, how- 
ever, are a different story: in all cases, Uranus acquires 
an eccentricity comparable to 0.1. Simulations reveal 
that the configurations where Saturn and Uranus are in 
a 3:2 MMR do not result in strong instabilities. This 
is because the system is given a chance to encounter 
high-order MMR's between Uranus and Neptune and 
spread out before crossing the 2:1 MMR between Sat- 
urn & Uranus. Furthermore, in order to increase the 
chances of Uranus/Saturn orbital crossing, it helps to 
start the two planets in the most compact stable reso- 
nance - namely the 4:3. From here, the degeneracy lies 
in whether Uranus and Neptune start out in a 3:2 or 4:3 
MMR. Out of 30 numerical simulations performed for 
each configuration, we only observed ice-giant/gas-giant 
scattering events in the evolutions of the system where 
Uranus and Neptune are in the 4:3 MMR. Particularly, 
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Fig. 8. — Dynamical evolution of the initial configuration where 
initially Jupiter & Saturn are in a 2:1 MMR, while both Saturn 
&: Uranus and Uranus & Neptune are in a 4:3 MMR's. In this 
model, the global instability is brought forth by Saturn & Uranus 
encountering a mutual 3:2 MMR. All else is as in fig.l. 

23% were successful, and in 57% of the successful in- 
tegrations, an ice giant exhibited encounters with both 
gas giants. Figure (8) shows the dynamical evolution of 
this configuration with time, while table (2) lists the final 
values of semi-major axes, eccentricities, and inclinations 
for the planets. 

Aside from Saturn's close encounter with an ice-giant, 
Jupiter's and Saturn's migration is dominated by scat- 
tering of planetesimals. As a result, equation (4) approx- 
imates their evolution well. Similarly to the previous sec- 
tion, when applied to Jupiter and Saturn simultaneously, 
this relation can be used to roughly estimate the starting 
semi-major axis of Jupiter and the initial rridisk- Setting 
Aa s = (2) 2 / 3 ay - (5/2) 2 / 3 a^ and Aaj = a f j - a l 3 with 
a j — 5.2 AU yields a total scattered mass of 49m® and 
cij = 5.5 AU. These values are in good agreement with 
numerical integrations. For instance, the evolution pre- 
sented in figure (1) started had the parameters: a l j = 5.5 
AU and m disk = 51m© . 

4. DISCUSSION 

The calculations presented in this work aim to place 
constraints on the early dynamical state of the solar sys- 
tem. We began by constructing a large array of multi- 
resonant systems. We then considered each one of these 
systems individually and tested them against two inter- 
related constraints: scattering of an ice-giant from a gas- 
giant and the outer ice giant undergoing a transient phase 
of high eccentricity. Both events are needed to reproduce 
the current dynamical architecture of the solar system. 

We showed numerically that three configurations in 
which Jupiter & Saturn start out in a 3:2 MMR follow 
evolutionary tracks which are compatible with our con- 
straints. The three systems are those in which (1) Saturn 
& Uranus are in a 3:2 MMR while Uranus & Neptune are 
in a 5:4 MMR, (2) Saturn & Uranus are in a 4:3 MMR 
while Uranus & Neptune are in a 3:2 MMR and (3) Sat- 
urn & Uranus are in a 4:3 MMR while Uranus & Nep- 
tune are also in a 4:3 MMR. Two of these configurations 
were previously thought to be unstable (Morbidelli et al. 
2007). 



We showed that the classic Nice model can be re- 
produced from multi-resonant initial conditions where 
Juipter & Saturn start out in a 5:3 MMR. The four 
particular primordial states that we found were those 
where initially, (1) Saturn & Uranus are in a 3:2 MMR 
while Uranus & Neptune are in a 5:4 MMR, (2) Saturn 
& Uranus are in a 3:2 MMR while Uranus & Neptune are 
in a 6:5 MMR, (3) Saturn & Uranus are in a 4:3 MMR 
while Uranus & Neptune are in a 3:2 MMR, and (4) Sat- 
urn & Uranus are initially in a 4:3 MMR while Uranus 
& Neptune are also in a 4:3 MMR. 

Finally, we used an analytical technique to rule out a 
large portion of the generated multi-resonant systems, in 
which Jupiter and Saturn are initially in the 2:1 MMR, 
based on an argument that the eccentricities generated 
by any Uranus/Neptune resonant encounters are too 
small. Simultaneously, we showed that the considered 
constraints can be satisfied by the dynamical evolution 
whose initial condition has Jupiter & Saturn locked in 
a 2:1 MMR, and the other pairs of the planets in 4:3 
MMR's. 

The calculations presented here are intended in part as 
a point of departure for future research. There is indeed 
a large array of unexplored issues. For instance, in the 
case of the system where Jupiter & Saturn start out in a 
2:1 MMR, it is not clear if a single scattering event alone 
is enough to correctly reproduce the secular dynamics of 
Jupiter and Saturn, or if a resonant encounter is also re- 
quired. In the initial condition where Jupiter and Saturn 
start out in a 5:3 MMR, dynamical stability may pose an 
issue (Morbidelli & Crida 2007). A further criterion of 
interest is LHB. While we do not attempt to time the 
onset of instability with LHB in these simulations, it is 
certainly fair to ask if the dynamical evolutions presented 
here are compatible with a long quiescent period preced- 
ing any resonant encounters. Two related issues imme- 
diately follow. First, will changing the placement of the 
inner boundary of the planetesimal swarm qualitatively 
change the process responsible for the onset of the insta- 
bility of the dynamical evolutions? Second, how will the 
effective viscosity that arises from the self-gravity of the 
disk affect our results? A large-scale numerical modeling 
effort will be instrumental in providing these answers. 

There is certainly room for broader study of the cur- 
rent setup of the problem as well. In this work, we have 
restricted ourselves to four-planet multi-resonant config- 
urations. Certainly, the idea of initially forming more 
than two ice giants is not unreasonable (Ford & Chiang 
2007, see however Levison & Morbidelli 2007). Although 
the results of Morbidelli et al. (2007) suggests that addi- 
tional planets in a compact multi-resonant system com- 
promise dynamical stability, more work is needed to ob- 
tain a good handle on this part of the problem. 

In conclusion, the determination of a small subset of 
initial conditions allows for a much more efficient survey 
of the parameter space. In this work, we have taken a 
step in this direction. We must keep in mind that the 
system at hand is highly chaotic, and must in the end be 
studied numerically. The resulting determinations are 
often probabilistic rather than conclusive, however the 
results are certainly bound to gain statistical weight as 
the number of completed simulations increases. Thus, 
while much progress is yet to be made, additional re- 
search carries great value since a solid understanding of 
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TABLE 2 

Orbital Elements at the End of Simulations 
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Orbital elements of solar system analogues, resulting from different initial conditions, at the end of the dynamical evolution simulations, 
presented in figures 1-8. Simulations where the ice giants switched places are labeled with a f . The disk mass used in each simulation, as 
well as the disk's power law index k are also given. 

initial conditions plays an unavoidably important role in We thank Hal Levison, Alessandro Morbidelli, Ramon 

further development of a comprehensive model for solar Brasser, Gregory Laughlin and Darin Ragozzine for use- 
system's formation. ful discussions. 
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TABLE 3 

Analytical Calculation of Planetesimal-Driven Migration 



Planet 


/ 


a (AU) (J:S 2:1) 


a (AU) (J:S 5:3) 


a (AU) (5:4 4:3 4:3) 


Jupiter 


0.35 


5.45 


5.69 


5.7 


Saturn 


0.2 


8.74 


7.97 


8.02 


Uranus 


0.08 


12.61 


9.84 


9.71 


Jupiter 


0.15 


18.36 


12.34 


11.77 



Inputs and results of the analytical calculations of planetesimal-driven migration (section 3.2). / is the effective scattering efficiency, 
inferred from numerical simulations. Although it varies from run to run, the variation is not too great. The second column, a (AU) (J:S 
2:1), lists the positions of the planets at the time of Jupiter/Saturn 2:1 MMR crossing, as inferred from the results of Levison et al. (2008). 
The next column, a (AU) (J:S 5:3), list the semi-major axes of the planets with Jupiter & Saturn nominally in the 5:3 MMR, traced 
back using equation (4) from the previous column. The last column, a (AU) (5:4 4:3 4:3), lists the semi-major axes of the 5:4 4:3 4:3 
multi-resonant configuration, assembled as discussed in section 2. Note the quantitative similarity between the two right-most columns. 
This leads us to believe that the 5:4 4:3 4:3 multi-resonant configuration, as well as other compact configurations from the same family are 
compatible with the classical Nice model. 



TABLE 4 

Analytical Estimates of Eccentricities After a Resonant Encounter 
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